Homogeneous cooling of rough, dissipative particles: Theory and simulations 
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We investigate freely cooling systems of rough spheres in 
two and three dimensions. Simulations using an event driven 
algorithm are compared with results of an approximate ki- 
netic theory, based on the assumption of a generalized homo- 
geneous cooling state. For short times t, translational and 
rotational energy are found to change linearly with t. For 
large times both energies decay like t~ 2 with a ratio indepen- 
dent of time, but not corresponding to equipartition. Good 
agreement is found between theory and simulations, as long 
as no clustering instability is observed. System parameters, 
i.e. density, particle size, and particle mass can be absorbed 
in a rescaled time, so that the decay of translational and rota- 
tional energy is solely determined by normal restitution and 
surface roughness. 

PACS: 51.10.+y, 46.10.+Z, 05.60.+W, 05.40.+j 



I. INTRODUCTION 

Collections of macroscopic, dissipative, and rough con- 
stituents have attracted a lot of interest recently, mainly 
in the context of granular media In the so-called 

grain-inertia regime, one is concerned with rapid granu- 
lar flow, which on one hand can be described by methods 
of statistical mechanics, analogous to the kinetic theory 
of dense gases and, on the other hand, can be well 

simulated with help of event-driven algorithms |p|,[lOj. In 
both approaches the dynamics of the system is assumed 
to be dominated by two-particle collisions which are mod- 
elled by their asymptotic states. A collision is character- 
ized by the velocities before and after the contact, and 
the contact is assumed to be instantaneous. In the sim- 
plest model, one describes inelastic collisions by normal 
restitution only. However, surface roughness is impor- 
tant for the dynamics of granular flow J8|— |To|] and allows 
for an exchange of translational and rotational energy. 

The rough sphere model was first introduced for elasti- 
cally colliding particles, assuming either perfectly rough 
or perfectly smooth particles |lj],[r2||. Subsequently, it 
has been generalised to intermediate values of the rough- 
ness to account for tangential restitution in inelastic col- 
lisions. Several groups have investigated the exchange of 
translational and rotational energy, using kinetic theory 
or numerical simulations One major result 

is that in contrast to conservative systems, energy is not 
equipartitioned between the degrees of freedom in dissi- 



pative systems [|4j— 11 0[| , even though the ratio of transla- 
tional and rotational energy approaches a constant, while 
both functions decay to zero in a freely cooling system. 
In this paper we study the full time evolution of trans- 
lational and rotational energy of a gas of rough, hard 
spheres. We compare numerical simulations to kinetic 
theory, which is based on a pseudo-Liouville-operator 
formalism and makes use of the assumption of a homo- 
geneous state. The computation covers the full range of 
times, from the initial linear change of energies with time, 
to the asymptotic state for large times, and also includ- 
ing the crossover between the two limiting regimes. 

In section [n] we introduce the microscopic interaction 
laws and the operators needed in section III to derive the 



solution for homogeneously cooling systems. The analyt- 
ical solution is compared to the simulations in section [Tv] 
and the results are summarized and discussed in |v|. 



II. MICROSCOPIC DYNAMICS 

We consider a system of N D-dimensional spheres 
interacting via a hard-core potential confined to a D- 
dimensional volume V = L D with linear size L and 
D = 2 or 3. The spheres have mass m, radius a, and 
moment of inertia /. The degrees of freedom are the po- 
sitions r M (i), the translational velocities v^t), and the 
angular velocities for each sphere, numbered by 

fi = 1,2,... ,N. Here, the spheres are rough with con- 
stant normal and also constant tangential restitution. 
The translational and angular velocities after collision 
(primed quantities) are determined by the velocities be- 
fore collision (unprimcd quantities) so that 



V n " — {V t + V r ) , 
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1 + r 



2q + 2 
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(V t + V r ) , 
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a(2q + 2) 

l + (3 
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[r x (v t + v r )] , and 
[r x (v t + v r )} . 
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The unit vector r = (r M — r^)/|r M — r„| points from 
particle v to \i, along the line connecting the centres of 
mass. We have introduced abbreviations for the normal 
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velocity v n = [(v^ — v v ) ■ r)] ■ r, for the tangential ve- 
locity due to translational motion v t — — v v — v n 
and for the tangential velocity due to rotational motion 
v r = —a^ft + ujv) x r. The constants r and (3 character- 
ize normal and tangential restitution, and q = I /(ma 2 ) 
depends on the mass distribution inside a grain. 

The time evolution of a dynamical variable A(t) — 
A(T;t) = ^({r^t),^),^^)};*) is (for positive 
times) determined by the pseudo-Liouville-operator C + 



A{t) = e^p(iC+t)A(0) 



(5) 



The pseudo-Liouville-operator ||, [l3],|l4| consists of two 
parts C+ = Cq + C + . The first, Co, describes the undis- 
turbed motion of single particles 



Co = -i^v^ ■ V r 



(6) 



and the second, C + = 5 X^^y (z 1 ^) ; describes hard- 
core collisions of all pairs of particles (p, v) with 

c+(pv) = 

i(v» v ■ r^Bi-v^ ■ r^Sdr^l - 2a)(fe+„ - 1) . (7) 

The operator acts on the particles p and v only, 
and instantaneously replaces the translational and an- 
gular momenta just before the collision, by the corre- 
sponding values just after. Q(~v^ u -r^ v ) is the Heaviside 
step-function which is non-zero for approaching particles 
only. The term v ^ v ■ r^ u gives the relative velocity of two 
colliding particles, reflecting that fast particles collide 
more frequently. Finally, we have introduced the nota- 
tion v^v = v ll — v v , r^y = r^-rv, and = r ^ / \r ^ . 



III. HOMOGENEOUS COOLING STATE 

The ensemble average of a dynamical variable A(T;t) 
is defined by 

(A) t = J dT p(T; 0)A(T; t) = f dT p(T; t)A(T; 0) (8) 
with the abbreviation for phase space integration 



dT = 



11^ 



«|-2o). (9) 



Here p(T;t) is the TV-particle phase space distribution 
function, i.e. p(T; t) dT is the probability at time t to 
find particle 1 at (n, t»i, Wi), particle 2 at (r2, V2, W2), 
etc. The time evolution of the iV-particle distribution 
p(T; t) = exp (— iC\t) p(T; 0) is governed by the adjoint 

C^_ of the time evolution operator C+. 

The quantities of interest are the translational and ro- 
tational energies per particle E tI = rn/ (2N)J2 fi v 2 l and 



E m t — 1 1 (2iV) J2fj, ^fn as we ^ as the total kinetic energy 
E = E tI + E Iot . It is impossible to calculate these expec- 
tation values exactly. To make some progress we resort to 
an approximation, known as homogeneous cooling state. 
One assumes that the iV-particle distribution function 
is spatially homogeneous and depends on time only via 
the average kinetic energy of the grains. It was shown 
in H that initially equipartitioned translational and ro- 
tational energy change with different rates, suggesting a 
generalised cooling state of the form 



PhcsIT;*) ~ exp 



-N 



Et, 



E IC 



T tr (t) T rot (t) 



(10) 



The .ZV-particle distribution does not depend on the spa- 
tial coordinates {r^} as a consequence of the assumption 
of homogeneity. It depends on time only via the average 
translational and rotational energy: T tv (t) = 2/D(E tI ) t 
and T Iot (t) = 2/ (2D - 3) (E Iot ) t . The factors D and 
2D — 3 account for the number of translational and rota- 
tional degrees of freedom in D dimensions respectively. 
To study the time evolution of the average translational 
and rotational energy we consider the corresponding time 
derivatives: 



d_ 

df 



T tr (t) = ^ dTp(T-t)iC + E tI and 
i T «*® = J dT p(T;t)iC+E mt 



2D -3 



(11) 



Then, we assume the generalized homogeneous cooling 
state of Eq. (^o|), replacing p(T;t) by PHCs(r;i) in Eqs. 
(|TT|), and arrive at 



^ T tr(*) = -jj (iC + E tI ) ws and 
-jjT mt (t) = on — (i£ + E TOt ) HCS 



.> jj ... *-'""ia > • (12) 

All that remains to be done is a high dimensional phase 
space integral, the details of which are delagated to ap- 
pendix A, where we present calculations for D = 2. After 
integration over phase space has been performed, we find 



(iC+E tI ) ws = -GAT*' 2 + GBT^ z T mt and 

(iC+E mt ) ws = GBT t 3 / 2 - GCT t \ /2 T mt (13) 

with the constants A, B, C, and G, whose values depend 
on space dimensionality D. 

A. The 2D-Case 



,1/2, 



by 



In two dimensions the constants in Eqs. (]ll|) are given 



N Fk 1 — r 2 77 

G = 8a v\ ~9 2a ' A = —f- + i {1 - " ' 

V \ m 4 2 



c=JLfi-2 



2? V a, 



(14) 
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We have introduced the abbreviation 
r\ = q(l + (3)/(2q + 2), and g(2a) denotes the pair cor- 
relation function at contact. We consider homogeneous 
discs with q = 1/2. 

In a system of elastically colliding particles [jl2j , G T^J 2 
is just the Enskog collision frequency, whose inverse de- 
termines the mean time between collisions. We use 
this mean collision frequency of the elastic system to 

1 /2 

rescale time according to r = GTj (0)t, and further- 
more introduce the dimcnsionless translational tempera- 
ture T = T tr (t)/T tr (0) and the dimensionless rotational 
temperature R — T TOt (t)/T tr (0). Expressed in scaled 
units Eqs. @ and (]l|) read in 2D: 



—T = -AT 3 ' 2 + BT^R, 
dr 

-^-R = 2BT 3/2 - 2CT 1/2 R . 
dr 



(15) 



To find the asymptotic value of the ratio of K = T/R we 
consider the differential equation 



dT _ -AK 
dR ~ 



B 



2BK - 2C 
which can be solved by a constant 



(16) 



K = (2C-A+ v/(2C - A) 2 + 8B 2 /(4B) (17) 



in agreement with Ref. p0|| . 

Of particular interest is the long time decay of trans- 
lational and rotational energy as a function of r and (3. 
We would expect the decay to be slowest, if energy is lost 
only due to normal restitution and not due to roughness. 
Hence, as a function of f3, translational and rotational 
energy should persist for the longest times for (3 = ±1. 
In between, i.e. — 1 < f3 < 1, we expect a faster decay 
because surface roughness is an additional mechanism for 
the disspation of energy. These expectations are born out 
by the following detailed discussion of Eqs.(15 17]). We 
assume that the ratio K = T/R has reached its asymp- 



totic value given in Eq. 
tute R = T/K into Eq. 



for some t > tq and substi- 
and obtain 



Ji _ _prp3/2^ 



(18) 



The resulting equation is of the same functional form as 
for homogeneous cooling of smooth spheres, except for 
the coefficient F = A — B/K, which contains all the 
dependence on (3 and r. Its solution is given by 



T : 



T(r ) 



[1 + T(to)V2(F/2)(t-to)]' 



(19) 



Asymptotically, for r 
decays like T ~ (-^r)" 
R « (£y^ T) 



— » oo, the translational energy 
2 and the rotational energy like 

In a double-logarithmic plot of T and 



R against r one should observe straight lines with slope 
—2 and axial sections (at r — 1) given by —2 In (F/2) for 
the translational energy, and by — 2 In (F^/K/2) for the 
rotational energy. For the same slope a larger axial sec- 
tion implies persistence for longer times, i.e. a 'slower' 
decay in time. In Fig. [I] we plot the prefactors F and 
F\fK against j3 for r = 0.99. 
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FIG. 1. Dependence of the asymptotic decay- prefactors F 
and FsfK on (3 for r = 0.99. 

As a function of (3, F is smallest for f3 — *■ ±1 cor- 
responding to the cases where no energy is lost due 
to friction. The decay of translational energy is then 
pushed out to the longest timescales. In Fig. [y, the max- 
imum of F is reached for /3f r lax ~ 0.17, corresponding 
to the 'fastest' decay of T. For (3 = — 1 we find that 
F\f~K — 0, because the rotational energy remains con- 
stant as a function of time. F\f~K reaches its maximum 
for (3™^ re -0.29 f /3 t ™ ax , so that the decay of rota- 
tional energy is 'fastest' for a value of (3, different from 
the one where the translational energy decays fastest. In 
the limit (3 — > 1 the axial sections for T and R have ap- 
proximately the same values. In this case the ratio K 
is close to unity, reflecting the rather effective exchange 
of rotational and translational energy for rough spheres 
and thus approximate equipartition. 



B. The 3D-Case 

In 3D we consider spheres with q = 2/5. The constants 
in Eq. ( [T3| ) are given by 

G = 8(2a) 2 -J-g(2a), A = — h ry(l - ry), 

V V m 4 



B 



(I 



q \ q 



(20) 



As in two dimensions we introduce a dimensionless time 
r = |GT t 1 / 2 (0)i with the factor 2/D = 2/3, which ac- 
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counts for the number of translational degrees of freedom. 
The time dependence of the dimensionless translational 
temperature T = T tr (t)/T tr (0) and of the dimensionless 
rotational temperature R — T rot (i) /T tr (0) follows from 



d_ 

d_ 

7h 



T = -AT 3/2 + BT 1/2 R 
R = BT 3/2 - CT 1/2 R . 



(21) 



The asymptotic value of the ratio of K = T/R is given 
in Refs. § and @: 

K = (c - A + y/(C - A) 2 + AB 2 \ /(2B) . (22) 



The long time limit can be discussed as in the 2D-case. 
We find a very similar result for the dependence of the 
asymptotic decay on r and (3. 

The full time dependent solution, obtained by numer- 
ical integration of Eqs. ( |l5| ) and (pi]), will be compared 
with simulations in the next section. 



IV. NUMERICAL EXPERIMENTS 

Since we are interested in the behavior of granular par- 
ticles cooling over several decades in time, we use an 
event driven (ED) method. In these simulations, the par- 
ticles follow an undisturbed translational motion until an 
event occurs. An event is either the collision of two par- 
ticles or the collision of one particle with a boundary of a 
cell (in the linked-cell structure) . The cells have no effect 
on the particle-motion here; they were solely introduced 
to accelerate the search for future collision partners in 
the algorithm (see also Appendix B). Using the velocities 
just before contact we compute the particles' velocities 
just after the contact following Eqs. (0-Q). In the ED 
method the contact duration is implicitly zero, matching 
well the corresponding assumption of instantaneous con- 
tacts used for the theory. We remark that ED algorithms 
run into problems when the time between events t ev gets 
too small. In dense systems with strong dissipation t ev 
may tend towards zero. As a consequence the so-called 
"inelastic collapse" can occur, i.e. the divergence of the 
number of events per unit time. The problem of the in- 
elastic collapse [|l5| |l6| , can be handled using restitution 
coefficients dependent on the time elapsed since the last 
event jl7],[l^]. For the contact that occurs at time Uj, 
one uses r = 1 and (3 = — 1 if at least one of the involved 
partners had a collision with another particle later than 
tij — t c . The time t c can be identified as a typical duration 
of a contact. The effect of t c on the simulation results 
is negligible for large r and small t c , what we checked 
in the simulations. Only for very small values r < 0.6 a 
cut-off time t c = 10~ 6 s was needed to avoid the inelas- 
tic collapse. For a more detailed discussion of the ED 
algorithm used, see Appendix B. 



Every simulation is first equilibrated with r = 1 and 
[3 = —1 until the velocity distribution is Maxwellian. 
Then the restitution coefficients are set to the selected 
values. In order to classify the systems used for the 
simulations we need to specify the number of parti- 
cles N and the volume fractions p2D = (N/V)na 2 and 
P3D = (N/V)(4/3)ira 3 in two and three dimensions, re- 
spectively. 



A. Results in 2D 

For short times we can solve Eqs. (|l^) analytically and 
get T = 1 — At and R = (2B)t. Hence our data can be 
effectively collapsed for short times by rescaling time ac- 
corrding to t — > At and rotational temperature accord- 
ing to R — > RA/(2B). We compare the theoretical re- 
sult [numerical solution of Eq. (|lq)] with simulations for 
various sets of parameters. The scaling collapses data 
for different particle number, particle radius and den- 
sity on the same master curves, since all these depen- 
dencies are hidden in the rescaled time r. The shape of 
the master curves depends only on the restitution coef- 
ficients r and (3. In order to check the dependence of 
the solution on system size we simulate a large system 
with N = 99856 particles and a small system with only 
N = 198. For the small system we consider two volume 
fractions p — 0.25 or p = 0.01, the latter correspond- 
ing to an extremely dilute system. For the large system 
we always use p — 0.25. For the larger volume fraction 
we get numerically 50.25 (2a) ~ 1.58 in perfect agreement 
with the 2D equivalent of the Carnahan-Stirling formula 



g(2a) 



1 - 7p/16 
(1-P) 2 



(23) 



from Ref . [|19| . Initially, the normalized energies are T = 
1 and R — for all data presented here. 

In Fig. H we plot T against rescaled time At for various 
simulations with fixed value of r = 0.99 and different par- 
ticle number N, volume fraction p, and tangential resti- 
tution (3. We observe that all simulations with the same 
(3 collapse, independent of the specific value of N or p. 
This indicates that the dimensionless units T and r are 
indeed the system's inherent "temperature" and "time" . 
The temperature T decays continuously from its initial 
value T — 1, following the function T = 1 — At for short 
times At < 0.1 and decaying like r -2 for long times. 
The monotonic dependence of T on (3 is due to our scal- 
ing of the horizontal axis, which includes factors of (3 via 
A. If we plot T against t we observe that the crossover 
between the initial linear decay and the asymptotic t~ 2 
behavior is nonmonotonic with (3: The time of crossover 
is largest for = — 1, then decreases and reaches a min- 
imum around /3", lax ~ 0.17 and then increases again for 
/3 — ► 1 [in agreement with Fig.[lJ. 

The simulations of the large system deviate from the 
theory and also from the simulations for the small system 
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at large times or small T. The deviation from the theo- 
retical curve is due to the density instability, i.e. clusters 
of particles form and the homogeneous cooling assump- 
tion fails. Additional simulations show that the deviation 
from theory occurs earlier for stronger dissipation due to 
either smaller r or \(3\. 




Ax 



FIG. 2. T as function of rescaled time At. Different 
symbols correspond to different simulations with TV = 99856, 
p = 0.25 (large), N = 198, p = 0.25 (small), and N = 198, 
p = 0.01 (dilute). The coefficient of restitution is r — 0.99 
and the value of j3 is given in brackets in the inset. The curves 
represent numerical solutions of Eqs. (|l5|). 

In Fig. H we plot RA/(2B) against normalized time 
At for the same set of parameters as in Fig. [| As 
expected from the solution for small At, we find that 
RA/(2B) is proportional to At for small times At < 0.1. 
In this regime, the initially 'cold' rotational degrees of 
freedom are activated due to the transfer of linear to 
angular momentum during collisions. After some equi- 
librium between translational and rotational tempera- 
ture is achieved, both degrees of freedom loose energy 
in the long time limit. Like the translational tempera- 
ture, also the rotational temperature is independent of 
N and p, only r and (3 determine the decay of rotational 
energy. As for the translational energy, the observed, al- 
most monotonic dependence of the crossover time for the 
rotational energy on (3 is due to our scaled units t — * At 
and R — > RAj (2B). If we plot R against r we find a sim- 
ilar nonmonotonic dependence of the crossover time on 
(3 as for the translational energy with however the mini- 
mum of rotational crossover time occuring at a different 
value of /?™ t ax w -0.29 [see Fig. |. 

The results of simulations are in very good agree- 
ment with the theoretical results, which are based on 
the generalized homogeneous cooling assumption for the 
A^-particle distribution function. Note that the agree- 
ment extends over many orders of magnitude in time 



and covers the initial increase of R, the crossover and 
the algebraic decay for long times. 




0.1 1 10 100 

Ax 



FIG. 3. RA/(2B) as function of rescaled time At. The 
simulations are the same as in Fig . |^. The curves represent 
the numerical solutions of Eqs. (113) . 

In Fig. H we present the ratio of T and RA/(2B) as a 
function of scaled time At for a selected set of param- 
eters which were used in the simulations shown in Fig. 

Data for two values (3 = 1.0 and (3 = —0.5 are shown 
in Fig. U(a). For the latter we compare systems at dif- 
ferent densities and observe that the dilute simulations 
are in perfect agreement with the theory, whereas the 
dense system with p = 0.25 deviates. The small system 
shows a slightly smaller ratio (2B/A) T/R, whereas for 
the large system, the ratio (2B/A) T/R is larger than ex- 
pected and eventually diverges for large At. This is again 
the regime where the homogeneous cooling assumption 
fails. Interestingly, the large simulation with N — 99856, 
p = 0.25, and (3 = 1.0 is in reasonable agreement with 
theory. From Fig. |^(b), we learn that the simulations of 
dilute systems are always in good agreement with the- 
ory. The ratio (2B/A) T/R is constant for large times, 
but there is no equipartition of energies in the transla- 
tional and rotational degrees of freedom. 
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FIG. 4. (2B/A)T/R as a function of rescaled time At 
for some simulations from Fig. ^. (a) The simulation of the 
large system with N = 99856, p = 0.25, and /9 = -0.5 (solid 
squares) shows no saturation of (2B/A) T/R which instead 
diverges for large At. (b) Here only simulations of the dilute 
system with N = 198, and p = 0.01 are compared to theory. 
All show saturation of the ratio (2B/A) T/R for large times. 



B. Results in 3D 

As in 2D, we can solve Eq. ( pl| ) analytically for short 
times, and get T = 1 — At and R = Bt. Hence we rescale 
time t —t At and rotational temperature R — > RA/ B. 

We simulated various systems characterized by volume 
fraction p and particle number N . One system with den- 
sity p — 0.087 and N — 1331 particles is denoted as 
medium; two other systems with lower density have the 
parameters p = 0.0021, N = 4096 dilute, and p = 0.0023, 
N = 68921 large. The abbreviation corresponds to the 
density, only for the large system one should read "di- 
lute and large". Finally, a system with higher density, 



i.e. p = 0.23 and N = 54872 dense is examined. To cal- 
culate the pair correlation function at contact, we use the 
Carnahan-Starling formula 



Apg(2a) = 



1 



P- 



P 



(l-p)3 



- 1 = 4p 



1-P/2 

(l-p)3 



(24) 



in 3D from Ref. j23| . Initially, the normalized energies 
are T = 1 and R = for all data presented here. 



Constant r — 0.99, variable f3 and p 

In Fig. |^ we plot T against normalized time At for 
r = 0.99 and various values of the tangential restitu- 
tion (3. We observe a very similar picture as in 2D and, 
again, reasonable agreement between theory and simula- 
tion over many orders of magnitude in time. For f3 < 0.5 
most of the dependence on f3 is taken into account by our 
scaling, r — ► At, so that the scaled data almost collapse 
for (3 < 0.5. 
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FIG. 5. T as function of rescaled time At in 3D. The 
symbols correspond to simulations with TV = 1331, p = 0.087 
(medium), r — 0.99, and different j3 as given in brackets in 
the inset. The curves represent numerical solutions of Eqs. 
( ^l[ ) with the three-dimensional constants from Eqs. (|2c|). 

In Fig. [fj] we compare simulations of different systems 
with the numerical solution of Eqs. pll Only the dense 
simulations deviate from the theoretical result. The scal- 
ing of time with At is rather successful for small and 
moderately dilute systems; in dense systems a deviation 
from the theory occurs for large At. 
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Constant j3 = —0.9, variable r and p 
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FIG. 6. T as function of rescaled time At in 3D from 
simulations with r — 0.99, and j3 — —0.9 or j3 = +0.9 as given 
in brackets. Different symbols correspond to simulations with 
N = 1331, p = 0.087 (medium), N = 4096, p = 0.0021 
(dilute), N = 68921, p = 0.0023 (large), and TV = 54872, 
p — 0.23 (dense). The curves are numerical solutions of Eqs. 
( pl[ ) with the three-dimensional constants from Eqs. (|2c|) . 

In Fig. @ we plot RA/B versus normalized time At for 
medium and dense systems. As in 2D, we find that RA/B 
increases proportional to At for small times (At < 0.1), 
reflecting the activation of initially 'cold' rotational de- 
grees of freedom due to collisions. This feature as well 
as the full time dependence is well reproduced by our 
theoretical analysis. 
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In Fig. H we present T, RA/B and the ratio of T and 
RA/B as a function of scaled time At for r — 0.6 and 
[3 = —0.9, where interesting structure is observed. The 
open symbols correspond to the medium system with 
p = 0.087, N — 1331. The data are in good agree- 
ment with the theoretical curves, whilst we obtain sub- 
stantial differences between theory and simulation in the 
case of the dense system, due to the density instability 
[not shown here]. For the medium system the loss of 
energy during collisions is predominantly due to normal 
restitution and only after the translational energy has 
decayed to a very small value (T < 10 -5 ) does one ob- 
serve the energy loss due to friction. The two regimes 
can be discussed analytically with help of Eqs. (^TJ). For 
intermediate times, when the translational energy is still 
appreciable, the equations can be simplified for almost 
smooth spheres, i.e. (/? w —1): 



— T = -AT 3/2 
dr 

4-R = -CT 1/2 R . 
or 



(25) 
(26) 



We have neglected terms of C((l + (3) 2 ) and approximate 
A » (1 - r 2 )/4 and C » 5(1 + /3)/U. The solution for T 
is that of smooth spheres, decaying like T(t) ~ (At/2)~ 2 
for large r. Substituting this result into the equation for 
R, we find R(t)/R{t q ) = (r/r )- Q with a = 2C/A. Here 
To is some intermediate timescale, larger than the time for 
the initial increase of R, but smaller than the timescale 
to reach the true asymptotic state. The above algebraic 
decay is shown in Fig. [s] as a straight dashed line with 
a w 0.396. Once the translational energy has decayed to 
a very small value as compared to the rotational energy, 
all terms in the differential equations for R and T are 
equally important. We then observe a crossover from a 
T~ a to a t~ 2 decay of the rotational energy. This true 
asymptotic state is characterized by a constant ratio Tf R 
and has been discussed above. The crossover between the 
two regimes shows up as a parallel shift for T [see Figs. || 
and [| , because the translational energy decays like r 
in both regimes, but with a different prefactor. 



FIG. 7. RA/B as function of rescaled time At. The data 
are selected situations from Fig . |H| The curves represent the 
numerical solutions of Eqs. (|2l]). 
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FIG. 8. T, RA/B, and (B/A)T/R for 3D simulations with 
j3 = —0.9 and r = 0.6 in the medium (open symbols) and in 
the dense (solid symbols) system from Fig. The thick lines 
give the solution of Eqs. (|2l|). 

When r is increased to a value close to unity, i.e. the 
elastic case, the intermediate time regime disappears, be- 
cause normal and tangential restitution are equally im- 
portant. This is demonstrated in Fig. ^ where we show 
T plotted against the normalized time At for medium 
density p = 0.087, (3 = —0.9 and different values of r, as 
given in the legend. 
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FIG. 9. T as function of rescaled time At. The density 
is p — 0.087, f3 — —0.9 and the restitution coefficient r is 
given in the insert. The curves represent again the numerical 
solutions of Eqs. (pl|). 



decay of smooth spheres [see Eq.fl25|)], which is indepen- 
dent of r, because we use scaled time At. In the true 
asymptotic regime, all curves have the same slope with 
however an axial section, which increases with decreasing 
r. Without scaling t with A the axial section decreases 
with decreasing r reflecting the more efficient dissipation 
of energy for smaller r. The agreement between theory 
and simulations is quite good for values of r as low as 
r = 0.6, and even for r = 0.2 only the crossover regime 
is not captured by theory. 



V. SUMMARY AND DISCUSSION 

Homogeneous cooling of colliding inelastic rough 
spheres has been investigated with numerical simulations 
and an approximate kinetic theory in two and three di- 
mensions as well. We have confirmed that surface rough- 
ness is an important characteristic of the grains, in so far 
as it determines the decay of translational energy, i.e. the 
rate of cooling. If energy loss due to small normal and 
tangential restitution are comparable, then one observes 
an initial linear change of translational and rotational en- 
ergy, followed by a crossover to the asymptotic regime, 
where both functions decay like t~ 2 . This regime is char- 
acterized by a constant ratio T/R, whose value depends 
on both, r and (3. The dependence on (3 is nonmonotonic, 
the ratio being smallest for (3 — ±1. This nonmonotonic 
dependence on [3 also holds for the crossover time, which 
is longest for f3 = ±1. If the coefficients of normal and 
tangential restitution are such that energy is lost mainly 
due to normal restitution, then we observe an interme- 
diate time regime, in between the initial linear change 
and the true asymptotic behaviour with constant ratio 
T/R. This intermediate regime is also characterized by 
an algebraic decay of translational and rotational energy: 
Translational energy decays like for smooth spheres (£~ 2 ), 
whereas rotational energy decays with an exponent that 
depends continuously on r and (3. 

The theoretical approach is based on the assumption of 
a generalised homogeneous cooling state: The ./V-particle 
distribution is assumed to depend on time only via the 
average energies of translation and rotation. Based on 
this assumption we computed the time decay of transla- 
tional energy and rotational energy without further ap- 
proximations. Good agreement with numerical simula- 
tions was found for a large range of timescales and param- 
eter sets, provided no density instability builds up. The 
initial linear change, the asymptotic t~ 2 behaviour, the 
crossover in between as well as the intermediate algebraic 
decay for almost smooth spheres - all these features are 
accurately reproduced by our theoretical ansatz. These 
findings certainly support the assumption of a homoge- 
neous cooling state and suggest to expand around the 
HCS state to study deviations from homogeneous cool- 
ing. 



In the intermediate time regime all curves follow the 
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APPENDIX A: DETAILED CALCULATION IN 2D 

In this appendix we explain, as an example, the main 
steps to calculate (iC + E tr ) KCS of Eq. ( |l3| ) in 2D. The 
expectation value is calculated with the TV-particle dis- 
tribution function, properly normalised 



PHCs(r;i) 



1 

W \2irT tI (t) 



N 



I 



exp 



-vl + 



2wT Ioi (t) 
I 



N/2 



N 

^ V2T tr (t)" 1 ' ' 2T mt {t)~ v 



(Al) 



The angular velocity is a scalar in two dimensions, but 
a vector in more than two dimensions. Free streaming 
does not change the energy, so we have to take into ac- 
count only the collison operator iC' + and we keep the 
abbreviation for T from Eq. (||) so that 

(iC' + E tl ) ws = 

If 1 N 

-2E dTpnc S (r;t)C + (ap)—J2mvl = 

-^E [dTp ws (T;t)C + (a(3)^(v 2 a + v 2 ) . (A2) 

The binary collision operator C+ (a(3) gives a contribution 
only, if either v = a or if v = (3. Next, we introduce two 
S- functions, 

~ ^ / dT J dR i dR i S ( R i - r ^(R2 - r p ) 



p ws (T;t)C + (af3)-(v 2 +v 2 a ) , (A3) 



which allows us to replace r a by Ri and rp by R 2 in 
C+ (a (3). Integration over all r M of the respective part of 
Eq. (A3) can then be performed and yields a factor 



/ II II 9 d^l -2a)5{R 1 -r a )5{R 2 -r p ) = 

V N - 2 g(R 12 ) (A4) 



^=1 



The pair correlation function g(R\ 2 ) depends on I-R12I = 
\Ri — R 2 \ only. Similarly integration over all velocities 



and angular velocities with index p, and a =/= p ft gives 
1 due to normalization. We can sum over N(N — 1) 
identical integrals and get 



(iC' + E tI ) ws 
(N — 1) ( m 



2V 2 \27TT tI (t)J 2vT mt (t) 



exp 



dujiduj2 dRi dR 2 dv\dv 2 



I 



2T tr (i)^^ 2T Iot (t)' 
g(r){v 12 ■ f)9 (-«ia • f) 5 (\r\ - 2a) AE tI . (A5) 

The loss of translational energy of two colliding particles 
is denoted by AE tY and given by 



AE tI = -[2 V (v-l)(v 2 i2 - (V12 ■ rf) - 

(l/2)(l-r 2 )(v 12 -r) 2 + 2 V 2 a 2 n 2 



(A6) 



Here we use the abbreviations r\ = ^q+a ' ^ = ^ L ° 1 
uo 2 ) I y/2, and Ri — R 2 = r = rr. To perform the remain- 
ing integrations we substitute 

tt = -^=(oj 1 + iv 2 ), uj = -j=(u)x- UJ2), (A7) 

V = -y=(v 1 +v 2 ), v = —( Vl -v 2 ), (A8) 



R\ — R 2 



R = R\ . 



(A9) 



The Jacobian determinant for the above transformation 
is 1. Integration over u), V and R can be done, which 
all give the value 1 due to normalisation. The resulting 
integral is 



(iC'+Etr) 



(N- l)m / 21 \ 1/2 



exp 



mv 



HCS 4nT tr (t)V \2nT rot (t) 

2 m 2 



dQdrdv 



2T tr (i) 2T rot (t) 



g(r)(vr)e (-v f) 8 (\r\ -2a) 



-[2 V (rj-l)(v 2 -(vr) 2 ) 



- (1/2)(1 - r 2 )(v ■ rf + 2i] 2 a 2 n 2 ] 

The integration over \r\ yields 2ag(2a). Choosing e.g. r 
to point along the x-eods, the integrals over linear and 
angular velocities can easily be done as moments of a 
Gaussian distribution. The result is independent of f, so 
that the integration over r gives 2tt. Finally we assume 
that 1 <C N, app roximate N w (N — 1) and obtain the 
result of Eq. @ 



APPENDIX B: ED ALGORITHM 

Simple ED algorithms update the whole system after 
each event, a method which is straightforward, but inef- 
ficient for large numbers of particles. In Ref. |2(J an ED 







algorithm was introduced which updates only those two 
particles which were involved in the last collision. For this 
a double buffering data structure is implemented, which 
contains the 'old' status and the 'new' status, each con- 
sisting of: time of event, position, velocities, and event- 
partner. When a collision occurs, the 'old' and 'new' sta- 
tus of the participating particles are exchanged. Thus, 
the former 'new' status becomes the actual 'old' one, 
while the former 'old' status becomes the 'new' one and is 
free for future calculations. This seemingly complicated 
exchange from information is carried out extremely sim- 
ple and fast by only exchanging the pointers to the 'new' 
and 'old' status respectively. The 'old' status of parti- 
cle i has to be kept in memory, in order to calculate the 
time of the next contact, ty, of particle i with any other 
object j which can change its status due to a collision 
with yet another particle. During the simulation this 
may be neccessary several times so that the predicted 
'new' status has to be modified. An object j is either a 
particle (j = 1, i — 1, i + 1, N) or a cell-wall (j 
= N + 1, ... ). The minimum of all is stored in the 
'new' status of particle i, together with the correspond- 
ing partner j. Depending on the implementation, also 
positions and velocities after the collision can be calcu- 
lated. This would be a waste of computer time, since 
before the time t%j , the predicted partners i and j might 
be involved into several collisions with other particles, so 
that we apply a delayed update scheme The min- 

imum times of event, i.e. the times which indicate the 
next event for a certain particle, are stored in an ordered 
heap tree, such that the next event is found at the top 
of the heap with computational effort of 0(1); changing 
the position of one particle in the tree from the top to 
a new position needs C(logiV) operations. The search 
for possible collision partners is accelerated by the use of 
a standard linked-cell data structure and consumes 0(1) 
of numerical ressources. In total, this results in numer- 
ical effort of C(iVlogiV) for N particles. For a detailed 
description of the algorithm see Ref. J2(J . Using all these 
algorithmic tricks, we are able to simulate up to 10 5 parti- 
cles within reasonable time on a small workstation (IBM 
P43/133MHz) pj. The particle number is limited by 
RAM-size (64MB) rather than CPU-power. 
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